library(stats)
suppressPackageStartupMessages(
library(fpp2)
)
library(forecast)
suppressPackageStartupMessages(
library(dynlm)
)
library(magrittr)
suppressPackageStartupMessages(
library(tseries)
)
L=read.csv('LaborForce.csv',header=T)
Labor = ts(L$LNU01327662,start=1992,frequency=12)
autoplot(Labor,main="US Labor Force Participation Rate",ylab="Percent",xlab="Year")

acf(Labor)

pacf(Labor)

t = time(Labor)
model1=dynlm(Labor~t)
plot1=ts(predict(model1),start=1992,frequency=12)
autoplot(plot1)+autolayer(Labor)

model2=dynlm(Labor~t+sin(2*pi*t/13)+cos(2*pi*t/13))
plot2=ts(predict(model2),start=1992,frequency=12)
autoplot(plot2) + autolayer(Labor)

model3=dynlm(Labor~t+I(t*cos(2*pi*t/15))+I(t*sin(2*pi*t/15)))
plot3=ts(predict(model3),start=1992,frequency=12)
autoplot(plot3) + autolayer(Labor)

model4=dynlm(Labor~t+sin(2*pi*t/13)+t*cos(2*pi*t/13))
plot4=ts(predict(model4),start=1992,frequency=12)
autoplot(plot4) + autolayer(Labor)

model5=dynlm(Labor~t+sin(2*pi*t/14)+t*cos(2*pi*t/14))
plot5=ts(predict(model5),start=1992,frequency=12)
autoplot(plot5) + autolayer(Labor)

model6=dynlm(Labor~t+cos(2*pi*t/12))
plot6=ts(predict(model6),start=1992,frequency=12)
autoplot(plot6) + autolayer(Labor)

model7=dynlm(Labor~t+t*cos(2*pi*t/15)+t*sin(2*pi*t/15))
plot7=ts(predict(model7),start=1992,frequency=12)
autoplot(plot7) + autolayer(Labor)

model8=dynlm(Labor~I(t*t)+t*cos(2*pi*t/13)+sin(2*pi*t/13))
plot8=ts(predict(model8),start=1992,frequency=12)
autoplot(plot8) + autolayer(Labor)

t_adj=t-1992
model28=nls(Labor~a+b*t_adj+
c*cos(2*pi*t_adj/p)+d*sin(2*pi*t_adj/p)+
e*t_adj*cos(2*pi*t_adj/p)+f*t_adj*sin(2*pi*t_adj/p)+
g*cos(4*pi*t_adj/p)+h*sin(4*pi*t_adj/p)+
i*t_adj*cos(4*pi*t_adj/p)+j*t_adj*sin(4*pi*t_adj/p),
start=data.frame(a=82,b=-0.29,p=14,
c=-0.0934,d=0.28,
e=0.031317,f=-0.035067,
g=0.01,h=0.01,
i=0.01,j=0.01))
plot28=ts(predict(model28),start=1992,frequency=12)
autoplot(plot28) + autolayer(Labor)

t_adj=t-1992
model29=nls(Labor~a+b*t_adj+
c*cos(2*pi*t_adj/p)+d*sin(2*pi*t_adj/p)+
e*t_adj*cos(2*pi*t_adj/p)+f*t_adj*sin(2*pi*t_adj/p)+
g*cos(4*pi*t_adj/p)+h*sin(4*pi*t_adj/p),
start=data.frame(a=82,b=-0.29,p=14,
c=-0.0934,d=0.28,
e=0.031317,f=-0.035067,
g=0.01,h=0.01))
plot29=ts(predict(model29),start=1992,frequency=12)
autoplot(plot29) + autolayer(Labor)

t_adj=t-1992
model30=nls(Labor~a+b*t_adj+
c*cos(2*pi*t_adj/p)+
e*t_adj*cos(2*pi*t_adj/p)+
g*cos(4*pi*t_adj/p),
start=data.frame(a=82,b=-0.29,p=14,
c=-0.0934,
e=0.031317,
g=0.01))
plot30=ts(predict(model30),start=1992,frequency=12)
autoplot(plot30) + autolayer(Labor)

model31=dynlm(Labor~t+
t*cos(2*pi*t/15.3)+t*sin(2*pi*t/15.3)+
cos(4*pi*t/15.3)+sin(4*pi*t/15.3))
plot31=ts(predict(model31),start=1992,frequency=12)
autoplot(plot31) + autolayer(Labor)

summary(model1)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.38458 -0.35754 0.00369 0.36341 1.41844
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 680.603408 7.270444 93.61 <2e-16 ***
## t -0.300534 0.003624 -82.92 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.537 on 334 degrees of freedom
## Multiple R-squared: 0.9537, Adjusted R-squared: 0.9535
## F-statistic: 6876 on 1 and 334 DF, p-value: < 2.2e-16
summary(model2)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/13) + cos(2 * pi *
## t/13))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.03689 -0.29964 0.00565 0.33782 1.10905
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 657.877547 6.356894 103.490 < 2e-16 ***
## t -0.289201 0.003169 -91.260 < 2e-16 ***
## sin(2 * pi * t/13) -0.258539 0.034085 -7.585 3.37e-13 ***
## cos(2 * pi * t/13) -0.372962 0.036977 -10.086 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4471 on 332 degrees of freedom
## Multiple R-squared: 0.9681, Adjusted R-squared: 0.9678
## F-statistic: 3356 on 3 and 332 DF, p-value: < 2.2e-16
summary(model3)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + I(t * cos(2 * pi * t/15)) + I(t *
## sin(2 * pi * t/15)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1393 -0.3780 -0.0018 0.3711 1.1479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.553e+02 7.144e+00 91.725 <2e-16 ***
## t -2.879e-01 3.562e-03 -80.843 <2e-16 ***
## I(t * cos(2 * pi * t/15)) 1.738e-04 1.967e-05 8.835 <2e-16 ***
## I(t * sin(2 * pi * t/15)) -2.015e-05 1.934e-05 -1.042 0.298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4841 on 332 degrees of freedom
## Multiple R-squared: 0.9626, Adjusted R-squared: 0.9622
## F-statistic: 2847 on 3 and 332 DF, p-value: < 2.2e-16
summary(model4)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/13) + t * cos(2 *
## pi * t/13))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08380 -0.29193 0.02731 0.29376 0.97821
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 664.145333 6.140604 108.156 < 2e-16 ***
## t -0.292303 0.003061 -95.507 < 2e-16 ***
## sin(2 * pi * t/13) -0.249624 0.032473 -7.687 1.73e-13 ***
## cos(2 * pi * t/13) 55.931373 9.444133 5.922 7.95e-09 ***
## t:cos(2 * pi * t/13) -0.028067 0.004708 -5.962 6.39e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4255 on 331 degrees of freedom
## Multiple R-squared: 0.9712, Adjusted R-squared: 0.9708
## F-statistic: 2787 on 4 and 331 DF, p-value: < 2.2e-16
summary(model5)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + sin(2 * pi * t/14) + t * cos(2 *
## pi * t/14))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.14586 -0.30482 0.02367 0.29644 0.90439
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 654.846505 6.154874 106.395 < 2e-16 ***
## t -0.287650 0.003068 -93.756 < 2e-16 ***
## sin(2 * pi * t/14) -0.107086 0.032726 -3.272 0.00118 **
## cos(2 * pi * t/14) 79.061920 8.956710 8.827 < 2e-16 ***
## t:cos(2 * pi * t/14) -0.039613 0.004466 -8.870 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4183 on 331 degrees of freedom
## Multiple R-squared: 0.9721, Adjusted R-squared: 0.9718
## F-statistic: 2888 on 4 and 331 DF, p-value: < 2.2e-16
summary(model6)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + cos(2 * pi * t/12))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.15791 -0.32242 0.02468 0.29195 1.14316
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 672.40224 6.09863 110.25 <2e-16 ***
## t -0.29643 0.00304 -97.50 <2e-16 ***
## cos(2 * pi * t/12) -0.42973 0.03537 -12.15 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4477 on 333 degrees of freedom
## Multiple R-squared: 0.9679, Adjusted R-squared: 0.9677
## F-statistic: 5020 on 2 and 333 DF, p-value: < 2.2e-16
summary(model7)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + t * cos(2 * pi * t/15) + t * sin(2 *
## pi * t/15))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.14894 -0.26691 0.05073 0.30602 0.90700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 650.602081 6.119585 106.315 < 2e-16 ***
## t -0.285535 0.003051 -93.597 < 2e-16 ***
## cos(2 * pi * t/15) -62.443236 8.039196 -7.767 1.02e-13 ***
## sin(2 * pi * t/15) 70.253341 8.698830 8.076 1.27e-14 ***
## t:cos(2 * pi * t/15) 0.031317 0.004007 7.816 7.39e-14 ***
## t:sin(2 * pi * t/15) -0.035067 0.004337 -8.085 1.19e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.412 on 330 degrees of freedom
## Multiple R-squared: 0.9731, Adjusted R-squared: 0.9726
## F-statistic: 2383 on 5 and 330 DF, p-value: < 2.2e-16
summary(model8)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ I(t * t) + t * cos(2 * pi * t/13) + sin(2 *
## pi * t/13))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.13848 -0.27107 0.05675 0.29772 0.91265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.797e+03 1.826e+03 -4.816 2.23e-06 ***
## I(t * t) -2.351e-03 4.538e-04 -5.180 3.87e-07 ***
## t 9.139e+00 1.821e+00 5.019 8.49e-07 ***
## cos(2 * pi * t/13) 3.448e+01 9.994e+00 3.451 0.000632 ***
## sin(2 * pi * t/13) -1.858e-01 3.362e-02 -5.526 6.67e-08 ***
## t:cos(2 * pi * t/13) -1.739e-02 4.981e-03 -3.491 0.000547 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4098 on 330 degrees of freedom
## Multiple R-squared: 0.9733, Adjusted R-squared: 0.9729
## F-statistic: 2409 on 5 and 330 DF, p-value: < 2.2e-16
summary(model28)
##
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + d * sin(2 *
## pi * t_adj/p) + e * t_adj * cos(2 * pi * t_adj/p) + f * t_adj *
## sin(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p) + h * sin(4 *
## pi * t_adj/p) + i * t_adj * cos(4 * pi * t_adj/p) + j * t_adj *
## sin(4 * pi * t_adj/p)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 81.721164 0.284669 287.074 < 2e-16 ***
## b -0.278969 0.018513 -15.069 < 2e-16 ***
## p 15.624579 2.453688 6.368 6.54e-10 ***
## c -0.509344 0.122191 -4.168 3.94e-05 ***
## d 0.297235 0.730436 0.407 0.6843
## e 0.057260 0.022125 2.588 0.0101 *
## f -0.003323 0.081301 -0.041 0.9674
## g 0.154910 0.241634 0.641 0.5219
## h 0.054556 0.084684 0.644 0.5199
## i -0.001288 0.014332 -0.090 0.9285
## j -0.002584 0.013551 -0.191 0.8489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.405 on 325 degrees of freedom
##
## Number of iterations to convergence: 17
## Achieved convergence tolerance: 6.825e-06
summary(model29)
##
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + d * sin(2 *
## pi * t_adj/p) + e * t_adj * cos(2 * pi * t_adj/p) + f * t_adj *
## sin(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p) + h * sin(4 *
## pi * t_adj/p)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 81.756853 0.063446 1288.595 < 2e-16 ***
## b -0.281288 0.004088 -68.812 < 2e-16 ***
## p 15.305461 0.525971 29.099 < 2e-16 ***
## c -0.488013 0.069050 -7.067 9.58e-12 ***
## d 0.205529 0.149410 1.376 0.169886
## e 0.053475 0.007593 7.042 1.12e-11 ***
## f 0.006954 0.016004 0.435 0.664199
## g 0.127955 0.038381 3.334 0.000955 ***
## h 0.045083 0.056791 0.794 0.427860
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4039 on 327 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 3.75e-06
summary(model30)
##
## Formula: Labor ~ a + b * t_adj + c * cos(2 * pi * t_adj/p) + e * t_adj *
## cos(2 * pi * t_adj/p) + g * cos(4 * pi * t_adj/p)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 81.825556 0.047067 1738.479 < 2e-16 ***
## b -0.284602 0.003076 -92.517 < 2e-16 ***
## p 16.165036 0.132575 121.931 < 2e-16 ***
## c -0.402314 0.066260 -6.072 3.48e-09 ***
## e 0.051564 0.004396 11.730 < 2e-16 ***
## g 0.099924 0.033334 2.998 0.00293 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4236 on 330 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 3.962e-06
summary(model31)
##
## Time series regression with "ts" data:
## Start = 1992(1), End = 2019(12)
##
## Call:
## dynlm(formula = Labor ~ t + t * cos(2 * pi * t/15.3) + t * sin(2 *
## pi * t/15.3) + cos(4 * pi * t/15.3) + sin(4 * pi * t/15.3))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.09435 -0.25024 0.04469 0.28124 0.95038
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.421e+02 6.201e+00 103.546 < 2e-16 ***
## t -2.813e-01 3.092e-03 -90.996 < 2e-16 ***
## cos(2 * pi * t/15.3) -2.235e+01 8.887e+00 -2.515 0.012373 *
## sin(2 * pi * t/15.3) -1.054e+02 8.534e+00 -12.356 < 2e-16 ***
## cos(4 * pi * t/15.3) -1.281e-01 3.309e-02 -3.871 0.000131 ***
## sin(4 * pi * t/15.3) 4.461e-02 3.320e-02 1.343 0.180046
## t:cos(2 * pi * t/15.3) 1.104e-02 4.428e-03 2.494 0.013121 *
## t:sin(2 * pi * t/15.3) 5.274e-02 4.257e-03 12.389 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4032 on 328 degrees of freedom
## Multiple R-squared: 0.9743, Adjusted R-squared: 0.9738
## F-statistic: 1780 on 7 and 328 DF, p-value: < 2.2e-16
AIC(model1,model2,model3,model4,model5,model6,model7,model8,
model28,model29,model30,model31)
BIC(model1,model2,model3,model4,model5,model6,model7,model8,
model28,model29,model30,model31)
checkresiduals(model1)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 266.34, df = 24, p-value < 2.2e-16
checkresiduals(model2)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 242.49, df = 24, p-value < 2.2e-16
checkresiduals(model3)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 251.85, df = 24, p-value < 2.2e-16
checkresiduals(model4)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 233.31, df = 24, p-value < 2.2e-16
checkresiduals(model5)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 230.37, df = 24, p-value < 2.2e-16
checkresiduals(model6)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 243.2, df = 24, p-value < 2.2e-16
checkresiduals(model7)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 229.23, df = 24, p-value < 2.2e-16
checkresiduals(model8)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 228.81, df = 24, p-value < 2.2e-16
checkresiduals(model28)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model29)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model30)
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

checkresiduals(model31)

##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 232.38, df = 24, p-value < 2.2e-16